Lindemann Criterion and the Anomalous Melting Curve of Sodium 



O 

o 

(N 

G 
3 



M. Martincz-Canaleifl and A. BergaraQ 

Materia Kondentsatuaren Fisika Saila, Zientzia eta Teknologia Fakultatea, 
Euskal Herriko Unibertsitatea, 644 Postakutxatila, 48080 Bilbo, Basque Country, Spain 
Donostia International Physics Center (DIPC), Paseo de Manuel Lardizabal, 20018, Donostia, Basque Country, Spain ana 
Centro Mixto CSIC-UPV/EHU, 1072 Posta kutxatila, E-20080 Donostia, Basque Country, Spain 

Recent reports of the melting curve of sodium at high pressure have shown that it has a very 
steep descent after a maximum of around 1000K at 31 GPa. This is not due to a phase transition. 
According to the Lindemann criterion, this behaviour should be apparent in the evolution of the 
Debye temperature with pressure. In this work, we have performed an "ab-initio" analysis of the 
behaviour of both the Debye temperature and the elastic constants up to 102 GPa, and find a clear 
trend at high pressure that should cause a noticeable effect on the melting curve. 
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At low pressures, alkali metals adopt high-symmetry 
simple structures and exhibit nearly-free-elcctron (NFE) 
behaviour. Thus, physics textbooks have always de- 
scribed group I elements as simple metals. Recent work, 
however, has shown that, at high pressure, there is a sig- 
nificant deviation from that simple NFE behaviourism, 
due basically to electronic shell collapse^. Apparently 
simple behaviour is not restricted to electronic properties 
or crystal structures. The melting curves of the alkalis 
reported by BridgmanA in 1926 were also relatively sim- 
ple, fitting well to the monotonically increasing Simon's 
equation: 
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Bowever, this simple law breaks at high pressure, when 
T m does not follow eq. (1). Further experiments per- 
formed on the heavy alkalis Rb and Cs&l found maxima 
in the melting curves, at 5 GPa in Rb and, in Cs, a double 
maximum near the transition from Cs I to Cs II. Sodium 
shows an even more blatant deviation from eq. (1). Re- 
cent diamond anvil cell (DAC) experiments performed by 
Gregoryanz et al4 showed the striking behaviour of Na: 
T m had a maximum of ~ 1000A" around 31GPa, and it 
then steeply decreased to room temperature at 100 GPa. 

The usual ab-initio approach to melting point calcula- 
tion is based on phase coexistence simulations via first- 
principles molecular dynamics (FPMD). Bowever, this is 
complicated and computationally extremely expensive. 
Simpler models have been used successfully in the past 
with the goal of describing the melting curve of several 
material o 9 ' 10 . Bowever, these materials do not deviate 
significantly from the behaviour described by Simon's 
Law. Can a simpler calculation, based on semiempiri- 
cal criteria, reproduce the Na melting curve? This work 
is oriented towards answering this question, given that 
simple criterions such as the Lindemann criterion or the 
Born stability criterion need a very definite behaviour to 
satisfy the experimental melting curve of Na. 



The melting temperature calculations presented here 
will be based on the Lindemann melting criterion^. This 
model, based on the harmonic approximation, predicts 
that melting will occur when the root mean square dis- 
placement reaches a certain value (generaly about 1 /8th) 
of the mean interatomic distance. Numerically, 
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Although reasonable at first sight, this model has sev- 
eral defects. The assumption of harmonical forces dis- 
cards the bond breakage occurring in melting, and it 
does not take into account the liquid phase. On the 
contrary, both the liquid and solid phases are related to 
the melting temperature and its evolution with pressure 
via the Clausius-Clapeyron equation. Nevertheless, given 
the success various forms of the criterion have had (see 
e.g.—), it is certainly interesting to test it in sodium and 
its characteristic melting curve. 

It is not very wise to take the meaning of the con- 
stant in equation [2] very seriously, so the Lindemann cri- 
terion will be used here as a single-parameter model. Be- 
ing interested in using the least possible external knowl- 
edge, the free constant C can be calculated from a single 
point. Then, working out the evolution of T m requires 
the knowledge of T m (0 GPa) and the evolution of Qr>. 
Another approach, more useful to compare equation 2 
with experiment, would be a least squares fitting of C 
to the experiment. Bowever, the authors are more in- 
terested in using the least amount of parametres than in 
checking how good the fit of equation [2] is. 

The Debye temperature has been evaluated from the 
phonon spectra at different pressures in two different 
ways. Once the phonon dispersion is known, one could 
obtain Ojj from the phonon density of states (DOS). Ap- 
plying the Debye approximation to the harmonic crystal 
yields the well known DOS: 
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where ko is chosen to obtain the correct number of states. 
This condition is fulfilled if k D 3 = 6ir 2 N/V, where N 
represents the numer of atoms per unit cell and V is the 
volume of the unit cell. The average speed of sound c is 
obtained by fitting a line to the low cu region in a g{to) 
versus tu 2 plot. The Debye temperature is then easily 
obtained from its definition, ksGlD = hckjj. This is cor- 
rect, since the approximations made in the Debye model 
are negligible in the low w, T region. 

It is also possible to evaluate 9 d from the elastic con- 
stants. In a cubic crystal, the low £ frequencies in the 
(£,0,0) branch and the elastic constats are related in the 
following manner: 

2 2 

Cix = P^, C i4 = p^ r (4) 
In a similar fashion for the (£,£,0) branch, 



P (GPa) q(a u ) Cu C12 C44 9c (K) 

0.00 7.738 9.79 7.47 6.72 195.0 147.7 

10.53 6.60 43.9 36.9 14.1 284.5 229.0 

20.93 6.20 73.8 62.9 17.0 308.7 253.2 

26.73 6.05 88.8 76.8 17.9 316.2 260.5 

29.91 5.98 97.2 84.6 18.2 318.6 262.8 

43.11 5.75 131.2 115.3 18.7 323.5 260.1 
52.08 5.63 153.2 135.0 18.4 323.5 255.9 
63.50 5.50 182.7 162.4 17.6 319.4 249.8 
63.50 6.94 179.9 170.3 17.8 297.0 245.9 

73.12 6.82 204.8 194.7 16.3 285.1 234.1 
84.91 6.70 233.5 222.1 14.1 269.7 217.5 
93.83 6.62 255.1 242.4 12.1 256.2 203.8 
102.4 6.55 275.7 263.3 9.9 236.3 174.9 



TABLE I: Evolution with pressure of the elastic constants (in 
GPa) and Qd of Na. The 6d on the left has been calculated 
via the elastic constants, while the rightmost one is calculated 
from the phonon DOS. 



Cn — C12 — 2p~p-, C44 — P~^t~ 

C n + C 12 + 2C 44 = 2p^ (5) 

where LOti refers to the transversal mode contained in 
the z = plane. Symmetry considerations make these 
the only independent elastic constants in cubic crystals. 
This allows to readily obtain the bulk and shear moduli, 

B = ^(C 11 +2C 12 ) 

G = l(3C^ + C 1 i-C 12 ) (6) 
5 

Both moduli determine the long wavelegnth speed of 
longitudinal and transversal waves in a solid, via the fol- 
lowing equations: 

[g Ib + Ig 

In the Debye model of specific heat, one characterizes 
the behaviour of the material with 6n, related to the 
vibrational properties of a solid via an average sound 
velocity. The two can be related through the following 
relation, 

where the average sound velocity v m can be defined as 
follows: 
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So obtaining the Debye temperature by any of these 
two methods only require knowledge of the phonon spec- 
trum. This has been done fully "ab-initio ", using density- 
functional perturbation theory (DFP T 12 i 13 i 14 ) as imple- 
mented in the Quantum-ESPRESSO codei^. It is im- 
portant to note that phonons have been calculated in a 
regular mesh of the Brillouin zone. Then the interatomic 
force constant matrix has been obtained by interpolating 
the results to a Fourier series. This method has been pre- 
ferred because it then allows ready calculation of phonon 
DOS and frequencies for an arbitrary q. Convergence 
has been checked comparing the frequencies for selected 
q, where direct calculations were available. 



III. MELTING CURVE OF SODIUM 

To study sodium, an LDAii norm-conserving pseu- 
dopotential has been considered satisfactory. Conver- 
gence to the desired accuracy has been achievd with a 
30 Ry cutoff. Similarly, the summations over the Bril- 
louin zone have been performed in a, at worst, 18 18 
18 Monkhorst-Packi^ grid. In some cases, up to 22 22 
22 grids have been used. Another important point is 
the phonon grid used to interpolate the interatomic force 
constant (IFC) matrix. With the goal of getting a qual- 
ity IFC matrix, a T-centered 10 10 10 phonon grid has 
been used for both bec and fee Na. The quality of the 
IFC matrix has been checked by comparing the frequen- 
cies obtained directly from DFTP calculations and the 
interpolated values at selected g-points. 

With the approximations and numerical parameters 
fixed as mentioned above, the phonon spectrum of Na 
has been calculated from ambient pressure to 103 GPa, 
on the onset of the fee to c/16 phase transition. The 
phases probed have been bee (from to 63 GPa) and fee 
(63 to 103 GPa). Table Q] shows the calculated evolution 
of the elastic constants with pressure. 
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FIG. 1: Evolution with pressure of the Na shear modulus 
(solid line), together with the contribution of C44 (dashed 
line) and Cn — C12 (dotted line). 



The first thing one notices is that C44 reaches a maxi- 
mum at around 43 GPa and, after the phase transition, 
it decreases dramatically. Although C — \{C\\ — C12) 
increases, its growth can not compensate the decrease of 
C44, which causes the decay of the shear modulus. As a 
consequence of this, Qu has a maximum in this pressure 
range, which, according to equation [2J should also result 
in a maximum in the melting curve of sodium. 

We now turn our attention to the evolution of the 
phonon DOS, represented in figure[2j There, one can ap- 
preciate that the high-frequency region of the spectrum, 
consisting basically on longitudinal modes, increases with 
pressure. On the other hand, the low frequency van-Hove 
singularity corresponding to transversal modes decreases 
monotonically. From equations 01 [5] and [6] it is possible to 
appreciate that the shear modulus of fee Na is collaps- 
ing, a similar behvaviour to that described by Kechini^. 
What the Lindemann criterion suggests is that this has 
very definite consequences on the evolution of T m with 
pressure. From equation^ on the limit vi » Vt, the av- 
erage sound velocity v m will tend to Vt- Substituting v t 
into equation [8] results in ©£> oc \fGp~ 1 l & . The resulting 
decrease of ©£> , together with the theromdynamic stabil- 
ity condition of < means, according to equation [5] 
that the melting temperature will decrease. 

With these considerations in mind, the behaviour of 
the calculated melting curve is not surprising. The melt- 
ing curve, calculated from the clastic constants, is shown 
in figure [H The first thing to notice is how a simple crite- 
rion is able to predict qualitatively the behaviour of the 
melting curve of Na. It also predicts that, at 102 GPa, 
Na should be a liquid at ambient temperature. However, 
this method also has its shortcomings. The most notice- 
able is that, even though it predicts that the melting 
curve will have a maximum, it severely underestimates 
the maximum at 600K. Even if one takes into account 
the large error margins presented by Gregoryanz et al& 
and the known shortcomings of LDA, this is still a large 
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FIG. 2: Evolution with pressure of the phonon DOS in fee 
Na. A detail of the low uj region is also depicted. 
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FIG. 3: Melting curve of Na. The shown line is a least squares 
fit to the Kechin equation^ 



difference. In any case, it is encouraging that the theo- 
retical maximum is located at 25 GPa, very close to the 
pressure of the experimental maximum, 31 GPa. This 
difference can even be caused by the slow convergence of 
phonon frequencies for small q. 

It should be noted, however, that the free parameter 
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in the Lindcmann model has been determined from a 
single point. Another option would have been to use 
least-squares fitting to the experimental data, or even to 
the maximum of the curve, which would have resulted 
in better fits. This has not been done here because the 
authors were interested in using the least amount of ex- 
ternal information. 



Conclusion 

In this work we present an "ab-initio" analysis of the 
evolution of the phonon DOS and the elastic constants of 
Na with pressure. It has been shown that for pressures 
higher than 40 GPa the Debye temperature O £> and the 
shear modulus start decreasing. This decrease is even 
steeper after the bcc to fee transition. Additionally, an 
estimate of the melting curve of Na has been made, with 
the only knowledge of the melting temperature of sodium 
at ambient pressure. 

Another interesting conclusion one can draw from the 
results is that, although naive, the Lindcmann criterion 
supplemented with ab-initio phonon DOS correctly pre- 
dicts the existence of a maximum in the melting curve of 
sodium, at a pressure remarkably close to the experimen- 
tal figure. The predictions based on the criterion fail to 
predict a high enough melting temperature, although this 
could be aggravated by the slow convergence of phonons 
for small q. One might be worried that it cannot predict 
the flat slope of the melting curve on the beginning of the 



fcc-liquid part shown in 8 -. In any case, there were no ex- 
perimental points reported in that region. Furthermore, 
more elaborate calculations presented in^l also support 
the idea that < in that region. 

Although using the ®d from either the phonon DOS 
or the elastic constants does not result in qualitatively 
different results, it seems clear from Fig. [2] that an ap- 
proach based on (u>) will probably fail in sodium, as it 
increases with pressure. 

Given the qualitative agreement found for Na, one 
might be tempted to do some qualitative predictions for 
lithium, for which technical difficulties have slowed ex- 
perimental advance. Preliminary calculations^ show no 
collapse of the shear modulus, neither a decrease of the 
Debye temperature up to 35 GPa. Taking into account 
the underestimate of the melting temperature for sodium 
and that at low pressures T m should have a subtle posi- 
tive sloped, it seems most likely that T m will rise slowly 
until the bcc to fee phase transition, where a Rb-likc 
maximum might appear. 
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